Network Constrained SPPA

Import packages

pacman::p_load(sp, sf, spNetwork, tmap,tidyverse,spatstat,raster,maptools,kableExtra,plotly,ggthemes,onemapsgapi,devtools,rgdal)

The downloaded binary packages are in
    /var/folders/hz/d8slc6jj09353ck_sskyq8ww0000gn/T//Rtmpgmf740/downloaded_packages

Data

Geospatial Data

Import Data

beijing_sf <- st_read("data/Geospatial/beijing.geojson") %>%
  st_transform(crs=4555)

Reading layer `beijing' from data source 
  `/Users/annatrw/annatrw/IS415MAC_GRP12/analysis/Network-Constrained/data/Geospatial/beijing.geojson' 
  using driver `GeoJSON'
Simple feature collection with 16 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 115.4172 ymin: 39.44247 xmax: 117.5071 ymax: 41.05861
Geodetic CRS:  WGS 84

Select Column for Analysis

beijing_sf <- select(beijing_sf, neighbourhood, geometry)
beijing_sf

Simple feature collection with 16 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 115.4172 ymin: 39.44247 xmax: 117.5071 ymax: 41.05861
Geodetic CRS:  New Beijing
First 10 features:
          neighbourhood                       geometry
1                东城区 MULTIPOLYGON (((116.4423 39...
2                西城区 MULTIPOLYGON (((116.3916 39...
3                昌平区 MULTIPOLYGON (((116.0427 40...
4       大兴区 / Daxing MULTIPOLYGON (((116.7347 39...
5                房山区 MULTIPOLYGON (((116.2466 39...
6      怀柔区 / Huairou MULTIPOLYGON (((116.279 40....
7  门头沟区 / Mentougou MULTIPOLYGON (((115.563 39....
8        密云县 / Miyun MULTIPOLYGON (((116.8826 40...
9       平谷区 / Pinggu MULTIPOLYGON (((117.3813 40...
10     延庆县 / Yanqing MULTIPOLYGON (((116.279 40....

Check CRS Again

st_crs(beijing_sf)

Coordinate Reference System:
  User input: EPSG:4555 
  wkt:
GEOGCRS["New Beijing",
    DATUM["New Beijing",
        ELLIPSOID["Krassowsky 1940",6378245,298.3,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Topographic mapping."],
        AREA["China - onshore."],
        BBOX[18.11,73.62,53.56,134.77]],
    ID["EPSG",4555]]

Check if there is invalid geometry

length(which(st_is_valid(beijing_sf ) == FALSE))

[1] 0

Check Missing Value

beijing_sf[rowSums(is.na(beijing_sf))!=0,]

Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  New Beijing
[1] neighbourhood geometry     
<0 rows> (or 0-length row.names)

Aspatial Data

beforecovid <- read_csv("data/aspatial/beijing_beforecovid.csv")

aftercovid <- read_csv("data/aspatial/beijing_aftercovid.csv")

Select Column for Analysis

beforecovid <- subset(beforecovid, select=c(id, street,neighbourhood_cleansed,latitude, longitude))

aftercovid <- subset(aftercovid, select=c(id,neighbourhood,latitude, longitude))

Check Missing Value

sum(is.na(beforecovid$latitude))

[1] 0

sum(is.na(aftercovid$latitude))

[1] 0

sum(is.na(beforecovid$longitude))

[1] 0

sum(is.na(aftercovid$longitude))

[1] 0

Visualize the listings

beforecovid_sf <- 

st_as_sf(beforecovid, 
                          # our coordinates are the longitude and latitude
                          coords = c("longitude", 
                                     "latitude"), 
                          # the geographical features are in longitude & latitude, in decimals
                          # as such, WGS84 is the most appropriate coordinates system
                          crs = 4326) %>%st_transform(crs=4555)
  #afterwards, we transform it to SVY21, our desired CRS

aftercovid_sf <- st_as_sf(aftercovid, 
                          # our coordinates are the longitude and latitude
                          coords = c("longitude", 
                                     "latitude"), 
                          # the geographical features are in longitude & latitude, in decimals
                          # as such, WGS84 is the most appropriate coordinates system
                          crs = 4326) %>% st_transform(crs=4555)

The listing before covid

tmap_mode("plot")+
  qtm(beijing_sf) +
  tm_shape(beforecovid_sf)+
  tm_dots()

We could notice that there are quite a lot of listings locate at the central area of Beijing

The listing after covid

tmap_mode("plot")+
  qtm(beijing_sf) +
  tm_shape(aftercovid_sf)+
  tm_dots()

After covid, the listings has reduced significantly

Network Analysis of Road

Import the Road Data of Beijing and name it as ‘network’

network <- st_read(dsn="data/Geospatial/Network/shape",
                   layer="roads")%>% 
  st_transform(crs=4555)%>% filter(type %in% c("residential", "footway","path"))

Reading layer `roads' from data source 
  `/Users/annatrw/annatrw/IS415MAC_GRP12/analysis/Network-Constrained/data/Geospatial/Network/shape' 
  using driver `ESRI Shapefile'
Simple feature collection with 98252 features and 7 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 116.08 ymin: 39.68 xmax: 116.77 ymax: 40.17999
Geodetic CRS:  WGS 84

network_sf <- select(network,osm_id,type, geometry)
network_sf

Simple feature collection with 40274 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 116.08 ymin: 39.68 xmax: 116.7699 ymax: 40.17999
Geodetic CRS:  New Beijing
First 10 features:
     osm_id        type                       geometry
1   4822342 residential LINESTRING (116.3972 39.993...
2   4964849 residential LINESTRING (116.4645 39.953...
3   5168733     footway LINESTRING (116.4561 39.942...
4   5168864 residential LINESTRING (116.4559 39.945...
5   5212775     footway LINESTRING (116.4561 39.945...
6   5217979     footway LINESTRING (116.4561 39.937...
7   9701854 residential LINESTRING (116.4584 39.952...
8   9779360 residential LINESTRING (116.401 39.9116...
9   9932203 residential LINESTRING (116.4611 39.932...
10 14347213 residential LINESTRING (116.4424 39.935...

Check the CRS

st_crs(network_sf)

Coordinate Reference System:
  User input: EPSG:4555 
  wkt:
GEOGCRS["New Beijing",
    DATUM["New Beijing",
        ELLIPSOID["Krassowsky 1940",6378245,298.3,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Topographic mapping."],
        AREA["China - onshore."],
        BBOX[18.11,73.62,53.56,134.77]],
    ID["EPSG",4555]]

Check invalid geometry and missing value

length(which(st_is_valid(network_sf) == FALSE))

[1] 0

network_sf[rowSums(is.na(network_sf))!=0,]

Simple feature collection with 0 features and 2 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  New Beijing
[1] osm_id   type     geometry
<0 rows> (or 0-length row.names)

Focus on Three Study Area

1. Shijingshan District(石景山区)

shijingshan_sf <- beijing_sf %>% filter(neighbourhood == "石景山区")

plot(shijingshan_sf)

Only select the road that within Shijingshan District

network_shijingshan <- subset(network_sf, lengths(st_intersects(network_sf, shijingshan_sf))!=0,)

Only Select the listings that within Shijingshan District

shijingshan_beforecovidlisting <- subset(beforecovid_sf, lengths(st_intersects(beforecovid_sf, shijingshan_sf))!=0,)
shijingshan_aftercovidlisting <- subset(aftercovid_sf, lengths(st_intersects(aftercovid_sf, shijingshan_sf))!=0,)

Visualize the Listings and the Network in Shijingshan Before Covid

tmap_mode("plot")
tm_shape(shijingshan_sf) +
  tm_polygons() +
tm_shape(network_shijingshan) +
  tm_lines() +
tm_shape(shijingshan_beforecovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

Visualize the Listings and the Network in Shijingshan After Covid

tmap_mode("plot")
tm_shape(shijingshan_sf) +
  tm_polygons() +
tm_shape(network_shijingshan) +
  tm_lines() +
tm_shape(shijingshan_aftercovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

Assign CRS

crs <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
network_shijingshan<- st_transform(network_shijingshan, crs)
shijingshan_beforecovidlisting<- st_transform(shijingshan_beforecovidlisting, crs)
shijingshan_aftercovidlisting<- st_transform(shijingshan_aftercovidlisting, crs)

Network Constrained KDE (NetKDE) Analysis

Preparing the lixels objects

Before computing NetKDE, the SpatialLines object need to be cut into lixels with a specified minimal distance

lixels <- lixelize_lines(network_shijingshan, 
                         700, 
                         mindist = 350)

Generating line centre points

Next, lines_center() of spNetwork will be used to generate a SpatialPointsDataFrame (i.e. samples) with line centre points as shown in the code chunk below.

samples <- lines_center(lixels)

Performing NetKDE for Beijing Before Covid

densities <- nkde(network_shijingshan, 
                  events = shijingshan_beforecovidlisting,
                  w = rep(1,nrow(shijingshan_beforecovidlisting)),
                  samples = samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE)

samples$density <- densities
lixels$density <- densities

samples$density <- samples$density*1000
lixels$density <- lixels$density*1000

tmap_mode('view')
tm_shape(lixels)+
  tm_lines(col="density")+
tm_shape(shijingshan_beforecovidlisting)+
  tm_dots()

Network Constrained G- and K-Function Analysis

Performing NetKDE for Beijing After Covid

densities1 <- nkde(network_shijingshan, 
                  events = shijingshan_aftercovidlisting,
                  w = rep(1,nrow(shijingshan_aftercovidlisting)),
                  samples = samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg =10,
                  sparse = TRUE,
                  verbose = FALSE)

samples$density_after <- densities1
lixels$density_after<- densities1
samples$density_after <- samples$density_after*1000
lixels$density_after <- lixels$density_after*1000

tmap_mode('view')
tm_shape(lixels)+
  tm_lines(col="density_after")+
tm_shape(shijingshan_aftercovidlisting)+
  tm_dots()

2. Haidian 海淀

haidian_sf <- beijing_sf %>% filter(neighbourhood == "海淀区")%>%
  st_transform(crs=4555)

plot(haidian_sf)

Only select the road that within Haidian District

network_haidian <- subset(network_sf, lengths(st_intersects(network_sf, haidian_sf  ))!=0,)%>%
  st_transform(crs=4555)

haidian_beforecovidlisting <- subset(beforecovid_sf, lengths(st_intersects(beforecovid_sf, haidian_sf))!=0,) %>%
  st_transform(crs=4555)
haidian_aftercovidlisting <- subset(aftercovid_sf, lengths(st_intersects(aftercovid_sf, haidian_sf))!=0,)%>%
  st_transform(crs=4555)

Visualize the Listings and the Network in Haidian Before Covid

tmap_mode("plot")
tm_shape(haidian_sf) +
  tm_polygons() +
tm_shape(network_haidian) +
  tm_lines() +
tm_shape(haidian_beforecovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

Visualize the Listings and the Network in Haidian After Covid

tmap_mode("plot")
tm_shape(haidian_sf) +
  tm_polygons() +
tm_shape(network_haidian) +
  tm_lines() +
tm_shape(haidian_aftercovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

crs <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
network_haidian<- st_transform(network_haidian, crs)
haidian_beforecovidlisting<- st_transform(haidian_beforecovidlisting, crs)
haidian_aftercovidlisting<- st_transform(haidian_aftercovidlisting, crs)

Network Constrained KDE (NetKDE) Analysis

lixels_haidian <- lixelize_lines(network_haidian, 
                         700, 
                         mindist = 350)

samples_haidian <- lines_center(lixels_haidian)

samples_haidian<- st_transform(samples_haidian, crs)

densities_haidian <- nkde(network_haidian, 
                  events = haidian_beforecovidlisting,
                  w = rep(1,nrow(haidian_beforecovidlisting)),
                  samples = samples_haidian,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE)

samples_haidian$density <- densities_haidian
lixels_haidian$density <- densities_haidian

samples_haidian$density <- samples_haidian$density*1000
lixels_haidian$density <- lixels_haidian$density*1000

The listings before covid

tmap_mode('view')
tm_shape(lixels_haidian)+
  tm_lines(col="density")+
tm_shape(haidian_beforecovidlisting)+
  tm_dots()

After Covid Visualization

densities_haidian <- nkde(network_haidian, 
                  events = haidian_aftercovidlisting,
                  w = rep(1,nrow(haidian_aftercovidlisting)),
                  samples = samples_haidian,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE)

samples_haidian$density_after <- densities_haidian
lixels_haidian$density_after <- densities_haidian
samples_haidian$density_after <- samples_haidian$density_after*1000
lixels_haidian$density_after <- lixels_haidian$density_after*1000

tmap_mode('view')
tm_shape(lixels_haidian)+
  tm_lines(col="density")+
tm_shape(haidian_aftercovidlisting)+
  tm_dots()

FENGTAI 丰台

fengtai_sf <- beijing_sf %>% filter(neighbourhood == "丰台区 / Fengtai")%>%
  st_transform(crs=4555)
plot(fengtai_sf)

network_fengtai <- subset(network_sf, lengths(st_intersects(network_sf, fengtai_sf  ))!=0,)%>%
  st_transform(crs=4555)
plot(network_fengtai)

Only select the road that within Fengtai District

fengtai_beforecovidlisting <- subset(beforecovid_sf, lengths(st_intersects(beforecovid_sf, fengtai_sf))!=0,) %>%
  st_transform(crs=4555)
fengtai_aftercovidlisting <- subset(aftercovid_sf, lengths(st_intersects(aftercovid_sf, fengtai_sf))!=0,)%>%
  st_transform(crs=4555)

Visualize the Listings and the Network in Haidian Before Covid

tmap_mode("plot")
tm_shape(fengtai_sf) +
  tm_polygons() +
tm_shape(network_fengtai) +
  tm_lines() +
tm_shape(fengtai_beforecovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

Visualize the Listings and the Network in Haidian Before Covid

tmap_mode("plot")
tm_shape(fengtai_sf) +
  tm_polygons() +
tm_shape(network_fengtai) +
  tm_lines() +
tm_shape(fengtai_aftercovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

Network Constrained KDE (NetKDE) Analysis

Preparing the lixels objects

crs <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
network_fengtai<- st_transform(network_fengtai, crs)
fengtai_beforecovidlisting<- st_transform(fengtai_beforecovidlisting, crs)
fengtai_aftercovidlisting<- st_transform(fengtai_aftercovidlisting, crs)

lixels_fengtai <- lixelize_lines(network_fengtai, 
                         700, 
                         mindist = 350)

samples_fengtai <- lines_center(lixels_fengtai)

densities_fengtai <- nkde(network_fengtai, 
                  events = fengtai_beforecovidlisting,
                  w = rep(1,nrow(fengtai_beforecovidlisting)),
                  samples = samples_fengtai,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE)

samples_fengtai$density <- densities_fengtai
lixels_fengtai$density <- densities_fengtai
samples_fengtai$density <- samples_fengtai$density*1000
lixels_fengtai$density <- lixels_fengtai$density*1000

tmap_mode('view')
tm_shape(lixels_fengtai)+
  tm_lines(col="density")+
tm_shape(fengtai_beforecovidlisting)+
  tm_dots()

After covid

densities_fengtai_after <- nkde(network_fengtai, 
                  events = fengtai_aftercovidlisting,
                  w = rep(1,nrow(fengtai_aftercovidlisting)),
                  samples = samples_fengtai,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE)

samples_fengtai$density_after <- densities_fengtai_after
lixels_fengtai$density_after <- densities_fengtai_after
samples_fengtai$density_after <- samples_fengtai$density_after*1000
lixels_fengtai$density_after <- lixels_fengtai$density_after*1000

tmap_mode('view')
tm_shape(lixels_fengtai)+
  tm_lines(col="density_after")+
tm_shape(fengtai_aftercovidlisting)+
  tm_dots()

The graphs reveal that road segments (darker orange color) are the part of network with relatively higher density of the listings. Those roads are located at/near the edge of the district The patterns are applicable to all of the three districts (Shijingshan, Fengtai, Haidian)

We could observe that the dense central area has been affected the most, the regions that remain unaffected are located at the edge of the district, showing that the lockdown measures has greatly impact the most of listings, especially those in the center of the Beijing where more tourists attraction are located, those areas are also the area with strict lockdown policy implementation

Insights Gained

The measures implemented to curb the spread of the pandemic had significant geographical effects. The closure or restricted access to public attractions, mostly located in central areas, greatly impacted the home-sharing and accommodation businesses like Airbnb. Risk management in the hospitality sector should be aware of changes in geographical spread of lodging facilities and consider how to strategise and optimise future locations of such facilities to better serve the needs of the market.